#include <cmath>
#include <cfloat>		// Needed for maximum floating point value

#include "Feature.h"

#include <iostream>
using namespace std;

//////////////////////////////////////////////////////////////////
// Feature Object												//
//////////////////////////////////////////////////////////////////

// Constructor for feature.
Feature::Feature(int img_w, int img_h)
{
	_img_w = img_w;
	_img_h = img_h;
	GenerateRandomFeature();
}

// This function generates a new random feature.
void Feature::GenerateRandomFeature()
{
	// Generate four random numbers in the correct range for the 
	// top left and bottom right.
	int x1 = 0, x2 = 0, y1 = 0, y2 = 0;
	while(x1 == x2 || y1 == y2)
	{
		x1 = rand() % _img_w;
		x2 = rand() % _img_w;
		y1 = rand() % _img_h;
		y2 = rand() % _img_h;
		if(x1 > x2)
		{
			_tl_x = x2;
			_br_x = x1;
		}
		else
		{
			_tl_x = x1;
			_br_x = x2;
		}
		if(y1 > y2)
		{
			_tl_y = y2;
			_br_y = y1;
		}
		else
		{
			_tl_y = y1;
			_br_y = y2;
		}
	}

	// Generate the feature type
	if((rand() % 2) == 0)
		_f = EDGE_RATIO;
	else
		_f = DOMINANT_EDGE;

	// Generate the orientations, make sure they're different!
	_o1 = rand() % ORIENTATION_BINS;
	while(_o1 == _o2)
		_o2 = rand() % ORIENTATION_BINS;

	_f_max = 0;
	_f_min = FLT_MAX;
}

// This function returns an array of feature values for a
// specified integral histogram set.
float* Feature::EvaluateFeature(IntegralHistSet *histset)
{
	cout << "Evaluating features.." << endl;
	// Allocate memory for results
	float* f_vals = new float[histset->_n];
	float* f_vals_p = f_vals;

	// Iterate through integral histograms in the set and 
	// calculate the feature values.
	IntegralHistNode* temp = histset->_head;
	// Check what kind of feature needs to be computed.
	if(_f == EDGE_RATIO)
		// Compute edge ratio feature
		while(temp)
		{
			// Compute first edge feature
			float h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o1];// Top left
			float h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o1];					// Top right
			float h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o1];					// Bottom left
			float h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o1];											// Bottom right
			float e1 = h1 - h2 - h3 + h4;
			// Compute second edge feature
			h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o2];	// Top left
			h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o2];						// Top right
			h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o2];						// Bottom left
			h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o2];												// Bottom right
			float e2 = h1 - h2 - h3 + h4;
			// Check whether the minimum and maximum values need 
			// to be updated. Note that an epsilon variable has 
			// been used to prevent division by zero. It is also 
			// important to note that some values of e1 and e2 will 
			// be negative when they are in fact 0 due to floating 
			// point inaccuracies. The epsilon value stops small 
			// negative values from occuring, and if it is not 
			// used, other measures must be taken to compensate
			// for this.
			if((e1 + EPSILON)/(e2 + EPSILON) > _f_max)
				_f_max = (e1 + EPSILON)/(e2 + EPSILON);
			if((e1 + EPSILON)/(e2 + EPSILON) < _f_min)
				_f_min = (e1 + EPSILON)/(e2 + EPSILON);
			// Determine ratio
			*f_vals_p++ = (e1 + EPSILON)/(e2 + EPSILON);
			temp = temp->next;
		}
	else
		// Compute dominant orientation feature
		while(temp)
		{
			float etotal = 0;
			float e1 = 0;
			// Compute the edge features in sequence
			for(int k=0; k<ORIENTATION_BINS; k++)
			{
				float h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+k];	// Top left
				float h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+k];						// Top right
				float h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+k];						// Bottom left
				float h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+k];												// Bottom right
				etotal = etotal + h1 - h2 - h3 + h4;
				if(k == _o1)
					e1 = h1 - h2 - h3 + h4;
			}
			// Check whether the minimum and maximum values need 
			// to be updated.
			if((e1 + EPSILON)/(etotal + EPSILON) > _f_max)
				_f_max = (e1 + EPSILON)/(etotal + EPSILON);
			if((e1 + EPSILON)/(etotal + EPSILON) < _f_min)
				_f_min = (e1 + EPSILON)/(etotal + EPSILON);
			// Determine value
			*f_vals_p++ = (e1 + EPSILON)/(etotal + EPSILON);
			temp = temp->next;
		}
	cout << "Done evaluating features" << endl;
	return f_vals;
}

// This function returns an array of feature values for a
// specified integral ghost set.
float* Feature::EvaluateFeature(GhostSet *g_set)
{
	cout << "Evaluating features.." << endl;
	// Allocate memory for results
	cout << "This set has " << g_set->_n << endl;
	float* f_vals = new float[g_set->_n];
	float* f_vals_p = f_vals;

	// Iterate through integral histograms in the set and 
	// calculate the feature values.
	IntegralHistNode* temp = g_set->_head;

	// Check what kind of feature needs to be computed.
	if(_f == EDGE_RATIO)
		// Compute edge ratio feature
		while(temp)
		{
			// Compute first edge feature
			float h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o1];// Top left
			float h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o1];					// Top right
			float h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o1];					// Bottom left
			float h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o1];											// Bottom right
			float e1 = h1 - h2 - h3 + h4;
			// Compute second edge feature
			h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o2];	// Top left
			h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o2];						// Top right
			h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+_o2];						// Bottom left
			h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+_o2];												// Bottom right
			float e2 = h1 - h2 - h3 + h4;
			// Check whether the minimum and maximum values need 
			// to be updated. Note that an epsilon variable has 
			// been used to prevent division by zero. It is also 
			// important to note that some values of e1 and e2 will 
			// be negative when they are in fact 0 due to floating 
			// point inaccuracies. The epsilon value stops small 
			// negative values from occuring, and if it is not 
			// used, other measures must be taken to compensate
			// for this.
			if((e1 + EPSILON)/(e2 + EPSILON) > _f_max)
				_f_max = (e1 + EPSILON)/(e2 + EPSILON);
			if((e1 + EPSILON)/(e2 + EPSILON) < _f_min)
				_f_min = (e1 + EPSILON)/(e2 + EPSILON);
			// Determine ratio
			*f_vals_p++ = (e1 + EPSILON)/(e2 + EPSILON);
			temp = temp->next;
		}
	else
		// Compute dominant orientation feature
		while(temp)
		{
			float etotal = 0;
			float e1 = 0;
			// Compute the edge features in sequence
			for(int k=0; k<ORIENTATION_BINS; k++)
			{
				float h1 = ((_tl_y == 0)||(_tl_x == 0)) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+k];	// Top left
				float h2 = (_tl_y == 0) ? 0 : temp->_hist->_histdata[((_tl_y-1)*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+k];						// Top right
				float h3 = (_tl_x == 0) ? 0 : temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+((_tl_x-1)*ORIENTATION_BINS)+k];						// Bottom left
				float h4 = temp->_hist->_histdata[(_br_y*temp->_hist->_w*ORIENTATION_BINS)+(_br_x*ORIENTATION_BINS)+k];												// Bottom right
				etotal = etotal + h1 - h2 - h3 + h4;
				if(k == _o1)
					e1 = h1 - h2 - h3 + h4;
			}
			// Check whether the minimum and maximum values need 
			// to be updated.
			if((e1 + EPSILON)/(etotal + EPSILON) > _f_max)
				_f_max = (e1 + EPSILON)/(etotal + EPSILON);
			if((e1 + EPSILON)/(etotal + EPSILON) < _f_min)
				_f_min = (e1 + EPSILON)/(etotal + EPSILON);
			// Determine value
			*f_vals_p++ = (e1 + EPSILON)/(etotal + EPSILON);
			temp = temp->next;
		}
	cout << "Done evaluating features" << endl;
	return f_vals;
}

// This function returns a random threshold. The threshold is
// chosen uniformly at random in range between the smallest and
// largest feature values.
float Feature::GenerateThreshold()
{
	float thresh = _f_min + ((_f_max-_f_min)*(float)rand())/(float)RAND_MAX;
	return thresh;
}

// This funtion prints out data describing the random feature.
// Useful for debugging.
void Feature::PrintData()
{
	cout << "The feature is a ";
	if(_f == DOMINANT_EDGE)
	{
		cout << "dominant orientation feature" << endl;
		cout << "Orientation " << _o1 << " is used" << endl;	
	}
	else
	{
		if(_f == EDGE_RATIO)
		{
			cout << "ratio of edges feature" << endl;
			cout << "Orientations " << _o1 << " and " << _o2 << " are used" << endl;
		}
	}
	cout << "The top left corner is (" 
		 << _tl_x << ", " << _tl_y 
		 << "), and the bottom right corner is ("
		 << _br_x << ", " << _br_y << ")" << endl;
}
